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Abstract. The collective behavior of cortical neurons is strongly affected by the presence of noise 
at the level of individual cells. In order to study these phenomena in large-scale assemblies of neurons, 
we consider networks of firing-rate neurons with linear intrinsic dynamics and nonlinear coupling, 
belonging to a few types of cell populations and receiving noisy currents. Asymptotic equations as 
the number of neurons tends to infinity (mean field equations) are rigorously derived based on a 
probabilistic approach. These equations are implicit on the probability distribution of the solutions 
which generally makes their direct analysis diflicult. However, in our case, the solutions are Gaussian, 
and their moments satisfy a closed system of nonlinear ordinary differential equations (ODEs), which 
arc much easier to study than the original stochastic network equations, and the statistics of the 
empirical process uniformly converge towards the solutions of these ODEs. Based on this description, 
we analytically and numerically study the influence of noise on the collective behaviors, and compare 
these asymptotic regimes to simulations of the network. We observe that the mean field equations 
provide an accurate description of the solutions of the network equations for network sizes as small as 
a few hundreds of neurons. In particular, we observe that the level of noise in the system qualitatively 
modifies its collective behavior, producing for instance synchronized oscillations of the whole network, 
desynchronization of oscillating regimes, and stabilization or destabilization of stationary solutions. 
These results shed a new light on the role of noise in shaping collective dynamics of neurons, and 
gives us clues for understanding similar phenomena observed in biological networks. 
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Introduction. The brain is composed of a very large number of neurons interact- 
ing in a complex nonlinear fashion and subject to noise. Because of these interactions, 
stimuli tend to produce coherent global responses, often with high reliability. At the 
scale of single neurons the presence of noise and nonlincarities often results in highly 
intricate behaviors. However, at larger scales, neurons form large ensembles that share 
the same input and are strongly connected, and at these scales, reliable responses to 
specific stimuli may arise. Such population assemblies (cortical columns or cortical 
areas) feature a very large number of neurons. Understanding the global behavior of 
these large-scale neural assemblies has been a fociis of many investigations in the past 
decades. One of the main interests of large-scale modeling is to characterize brain 
function at the scale most non-invasive imaging techniques operate. Its relevance is 
also connected to the fact that anatomical data recorded in the cortex reveal its colum- 
nar organization at scales ranging from about 50/xm to Irnrn. These so-called cortical 
columns contain of the order of hundreds to hundreds of thousands of neurons, and 
are generally dedicated to specific functions. For example, in the visual cortex VI, 
they respond to preferred orientations in visual stimuli. These anatomical and func- 
tional organizations point towards the fact that a significant amount of information 
processing does not occur at the scale of individual neurons but rather corresponds to 
a mesoscopic signal arising from the collective dynamics of many interacting neurons. 

In the computational ncuroscience community, this problem has been mainly 
studied using approaches based on the statistical physics literature. Most models 
describing the emergent behavior arising from the interaction of neurons in large-scale 
networks have relied on continuum limits since the seminal works of Wilson and Cowan 
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and Amari [3l |4] [53l |54] . Such models represent the activity of the network through 
a global variable, the population-averaged firing rate, which is generally assumed to 
be deterministic. Many analytical properties and numerical results have been derived 
from these equations and related to cortical phenomena, for instance in the case of 
the problem of spatio-temporal pattern formation in spatially extended models (see 
e.g. [T7J [211 ESI E] ) ■ This approach implicitly makes the assumption that the effect 
of noise vanishes in large populations. 

However, increasingly many researchers now believe that the different intrinsic 
or extrinsic noise sources participate in the processing of information. Rather than 
having a pure disturbing effect there is the interesting possibility that noise conveys 
information and that this can be an important principle of brain function [^21. In 
order to study the effect of the stochastic nature of the firing in large networks, many 
authors strived to introduce randomness in a tractable form. A number of compu- 
tational studies that successfully addressed the case of sparsely connected networks 
of integrate-and-fire neurons are based on the analysis of large assemblies that fire in 
an asynchronous regime [TJ [51 113j . Because of the assumption of sparse connectivity, 
correlations of the synaptic inputs can be neglected for large networks. The resulting 
asynchronous irregular states resemble the discharge activity recorded in the cerebral 
cortex of awake animals [19 . 

Other models have been introduced to account for the presence of noise in neu- 
ronal networks, such as the population density method and related approaches |15j . 
allowing efficient simulation of large neuronal populations. In order to analyze the 
collective dynamics, most population density-based approaches involve expansions in 
terms of the moments of the resulting random variables, and the moment hierarchy 
needs to be truncated in order to get a closed set of equations, which can raise a 
number of technical issues (see e.g. [55]). 

Yet other models of the activity of large networks are based on the definition of 
a Markov chain governing the firing dynamics of the neurons in the network, where 
the transition probability satisfies a differential equation called the master equation. 
Seminal works of the application of such modeling for neuroscience date back to the 
early 90s and have been recently developed by several authors ^40', 123] . Most of these 
approaches are proved correct in some parameter regions using statistical physics tools 
such as path integrals [13] and Van-Kampen expansions [H]. They motivated a num- 
ber of interesting studies of quasicycles |10J and power-law distribution of avalanche 
phenomena 7J. In many cases the authors consider one-step Markov chains, imply- 
ing that at each update of the chain, only one neuron in the whole network either 
fires or stops firing, which raises biological plausibility issues. Moreover, analytical 
approaches mainly address the dynamics of a finite number of moments of the firing 
activity, which can also raise such issues as the well-posedness [33] and the adequacy 
of these systems of equations with the original Markovian model [37] . 

In the present study, we apply a probabilistic method to derive the limit behavior 
resulting from the interaction of an infinite number of firing-rate neurons nonlinearly 
interconnected. This approach differs from other works in the literature in that it 
relies on a description of the microscopic dynamics (neurons in the network), and 
does not make the assumption of a sparse connectivity. Our model takes into account 
the fact that cortical columns feature different populations. The approach consists in 
deriving the limit equations as the total number of neurons tends to infinity, based 
on results obtained in the field of large-scale systems of interacting particles. This 
problem has been chiefly studied for solving statistical physics questions, and has been 
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a very active field of research in mathematics during the last decades [Ml [22l ESI HI] . 
In general, the equations obtained by such rigorous approaches are extremely hard 
to analyze. They can be either seen as implicit equations in the set of stochastic 
processes, or as non-local partial differential equations on the probability distribution 
through the related Fokker-Planck equations. But in both cases, understanding the 
dynamics of these equations is very challenging, even for basic properties such as the 
existence and uniqueness of stationary solutions and a priori estimates |30| . It appears 
even more difficult to understand qualitatively the effects of noise on the solutions and 
to interpret them in terms of the underlying biological processes. 

In the present article we aim at answering some of these questions. In the case 
we address, the problem is rigorously reducible to the analysis of a set of ordinary 
differential equations. This is because the solution of the mean field equations is a 
Gaussian process. It is therefore completely determined by its first two moments which 
we prove to be the solutions of ordinary differential equations. This allows us to go 
much deeper into the analysis of the dynamical effects of the parameters, in particular 
those related to the noise, and to understand their influence on the solutions. The 
analysis of this Gaussian process also provides a rich amount of information about 
the non-Gaussian solution of the network when its size is large enough. 

The paper is organized as follows. In the first section we deal with the modeling, 
the derivation of the mean field equations and of the related system of ordinary 
differential equations. We then turn in section [2] to the analysis of the solutions of 
these equations and the influence of noise. We show in details how noise strongly 
determines the activity of the cortical assembly. We then return to the problem 
of understanding the behavior of finite-size (albeit large) networks in section [s] and 
compare their behavior with those of the solutions of the mean field equations (infinite- 
size network). The analysis of the network behaviors in the different regimes of the 
mean field equations provides an interpretation of the individual behaviors responsible 
for collective reliable responses. This has a number of consequences that are developed 
in section m 

1. Model and mean field equations. In all the article, we work in a complete 
probability space (fi, J^, P) assumed to satisfy the usual conditions. 

We are interested in the large scale behavior arising from the nonlinear coupling 
of a large number N of stochastic diffusion processes representing the membrane 
potential of neurons in the framework of rate models (see e.g. [THl HH])- Hence the 
variable characterizing the neuron state is its firing rate, that exponentially relaxes to 
zero when it receives no input, and that integrates both external input and the current 
generated by its neighbortj^] The network is composed of P neural populations that 
differ by their intrinsic dynamics, the input they receive and the way they interact 
with the other neurons. Each population a £ {1, . . . , P} is composed of neurons, 
and we assume that the ratio N^/N converges to a constant 6a in ]0, 1[ when the total 
number of neurons N becomes arbitrarily large. We define the population function 
p that maps the index i e {1, . . . N} of any neuron to the index a of the population 
neuron i belongs to: p{i) = a. 

For any neuron i in population a, the membrane potential V/ has a linear intrinsic 
dynamics with a time constant r^. The membrane potential of each neuron returns 
to zero exponentially if it receives no input. The neuron i in population a receives 



'Alternatively, our model can be seen as a hypercolumn model, each diffusion process character- 
izing the activity of a whole cortical column modeled by Wilson and Cowan equations. 



4 



J. TOUBOUL, G. HERMANN & O. FAUGERAS 



an external current, which is the sum of a deterministic part Ia(t) and a stochas- 
tic additive noise driven by N independent adapted Brownian motions (-B')i=i...Ar 
modulated by the diffusion coefficients Xa{t). This additive noise term accounts for 
different biological phenomena [2], such as intrinsically noisy external input, channel 
noise |52j produced by the random opening and closing of ion channels, thermal noise 
and thermodynamic noise related to the discrete nature of charge carriers. 

Neurons also interact through their firing rates, classically modeled as a sigmoidal 
transform of their membrane potential. These sigmoidal functions only depend on the 
population a the neuron i belongs to: Sa{V'^). The functions Sa : M i— M are assumed 
to be smooth (Lipchitz continuous), increasing functions that tend to at — oo and to 
1 at oo. The firing rate of the presynaptic neuron j, multiplied by the synaptic weight 
Jij, is an input current to the postsynaptic neuron i. We classically assume that the 
synaptic weight Jij is equal to Jp(i)p{j) / Np(^j) . In practice this synaptic weight ran- 
domly varies depending on the local properties of the environment. Models including 
this type of randomness are not covered in this paper. The scaling assumption is 
necessary to ensure that the total input to a neuron does not depend on the network 
size. 

The network behavior is therefore governed by the following set of stochastic 
differential equations: 



dV\t) 



E 

r-p{j)= 



Sp{V'{t)) dt + K{t)dBl (1.1) 



These equations represent a set of interacting diffusion processes. Such processes 
have been studied for instance by McKean, Tanaka and Sznitman among others j36l 

Uni m] ■ This case is simpler than the different cases treated in where the 
intrinsic dynamics of each individual diffusion is nonlinear. 

In such models, the propagation of chaos property applies. This means that, 
provided the initial condition is drawn independently from the same law for all neurons 
of each population, in the limit N — > oo, all neurons behave independently, and have 
the same law which is given by an implicit equation on the law of the limiting process. 
This convergence and the limit equations are the subject of the following theorem: 

Theorem 1.1. Let T > a fixed time. Under the previous assumptions, we 
have: 

(i). The process for i in population a, solution of equation (1.1), converges in 



law towards the process V" solution of the mean field implicit equation: 



dV^{t) = 



1 P 

K(<) + I^{t) + ■^-'P^ [Sp{V[3{t))\ 



dt + K{t)dB'^{t) 

(1.2) 



as a process for t G [0, T], in the sense that there exists (V^/)t>o distributed as 
(V'(")t>o such that 



E 



sup \v:-v:\ 

0<t<T 



< 



C{T) 



N 



where C(-) is a function of time depending on the parameters of the system. 
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As a random variable, it converges uniformly in time in the sense that: 

sup E[\v;-v;\] < ^ 

o<t<T vN 



where C does not depend on time. In equations \1.2\ , the processes (i?"(i))Q=i 
are independent Brownian motions. 



(ii). Equation (1.2 I has a unique (trajectorial and in law) solution which is square 
integrable. 

(Hi). The propagation of chaos applies, i.e. provided that the initial conditions 
of all neurons are independent the law of {V'^'-{t), . . . ,V'^''{t),t < T) for 
any fixed k > 2 and {ii, . . . ,ik) E N'^, converges towards i^p(i-^) (E) . ■ . <E) Vp{i^) 



when N — > oo, where we denoted the law of the solution of equation (1.2 I 
corresponding to population a. This means that {y'''^(t), . . . , V'^''(t)) become 
independent processes. 

The proof of this theorem essentiaUy uses results from the works of Tanaka and 
Sznitman, summarized in j44j . A distinction with these classical results is that the 
network is not totally homogeneous but composed of distinct neural populations. 
Thanks to the assumption that the proportion of neurons in each population is non- 
trivial {Na/N —J' 6a g]0, 1[), the propagation of chaos occurs simultaneously in each 
population yielding our equations, as it is shown in a wider setting in |49j . The main 
deep theoretical distinction is that the theorem claims a uniform convergence in time: 
most of the results proved in the kinetic theory domain show propagation of chaos 
properties and convergence results only for a finite time, and convergence estimates 
diverges as time increases. Uniform propagation of chaos is an important property as 
commented in [16;, and particularly in our case as we will further comment. Methods 
to prove uniformity are generally involved (see e.g. |37) where uniformity is obtained 
for certain models using a dual approach based on the analysis of generator operators). 
Due to the linearity of the intrinsic dynamics, we provide here an elementary proof 
of this property in our particular system. 

Proof. The existence and uniqueness of solutions can be performed in a classical 



fashion using Picard iterations of an integral form of equation ( 1.2 ) and a contraction 
argument. The proof of the convergence towards this law, and of the propagation 
of chaos can be performed using Sznitman's powerful coupling method (see e.g. |44] . 
method which dates back from Dobrushin [H]), that consists in defining different 



independent processes V"^ solution of equation (1.2) driven by the same Brownian 
motion {Bl)t as involved in the network equation (1.1 1, and with the same initial 
condition V^{0) as neuron i in the network. Once these processes have been defined, 
it is easy to see that for a neuron i belonging to population a: 



0=1 "'o '5 r-pU)=P ^ 



ds 



^The initial conditions are said to be chaotic. 
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We have, denoting by t the maximal value of (r^, /3 = 1 . . . P): 



\vi-v:\<K 



.-its)/ 



^ max \Vl -Vl\ds 



K' 



NJo ^ 



3 = 1. ..-N 

t N 
{t-s)/r, 



(1.3) 



where K — maxa'^p \Jap\L. L is the largest Lipschitz constant of the sigmoids 
{Si3, l3 — I . . . P), and K' = maxo, \Jap \ max^ N/Np, quantity upperbounded, for 
N sufficiently large, by maxa \Jai3 \ max^ 2/5^. 

Since the righthand side of (1.3) does not depend on the index i, taking the 
maximum with respect to i and the expected value of both sides of (1.3 1, we obtain 

E 



max \v;-v; 

1=1. ..N 



< K / e-(*-")/^E 



max |V?"V7| 

j = l...N 



ds 



K'E 



max 

a = l P 



ho' 



N 



(1.4) 



Since the random variables Spi^j){Vj) — E [<5'p(j) (V7)] are independent and centered, 
and using the fact that the sigmoids Sp take their values in the interval [0, 1], we have: 



E 



max 



N 



N 



.-its)/ 



1 ^ r / /■* 

— ^E max(^y^ e-(*-)/-(5,(,)(y/)-E[5,(,)(V'/)])d- 



N 



< 



-y 



7V2 



N 



By Cauchy-Schwartz inequality, we can upperbound the second term of the right- 
hand side of inequahty (1.4 1 hy t/VN. Therefore, defining M( = E [max^ |V/ — F/j], 
we have: 

Mt<K f e"(*"'')/"'M^ds + X'^ 
Jo Vl 

implying, using Gronwall's lemma, 

Mt < 



N 

This inequality readily yields the almost sure convergence of V/ towards as N goes 
to infinity, uniformly in time, and hence convergence in law of V/ towards T^". 

The almost sure convergence of {Vl)t^[o,T] (considered as a process) towards 
iVi)te\o. T] ca n be proved in a similar fashion. Indeed, upperbounding the exponential 
term in (1.3) by 1 and taking the supremum, it is easy to see that: 



E 



sup Taax\V;-V;\ 

0<t<T«=l---^ 



< K I E 



sup niax - F/| 



K'T 



Noise effects in mean field dynamics 



7 



using the fact that: 



E 



N 



max sup , 

" tG[0,T] \ ^ Jo 



T f \„ 
- N^Jo 



< 



T 
N 



N 



^(^p(,)(?/)-E[5p(,)(Vy)]) 



^pO)(V;^)-E[^pO)(V'/)] 



using the independence of the and Cauchy- Schwartz inequahty. This last estimate 
readily implies, using Gronwall's inequality: 



E 



sup max \Vf - VI 

0<t<T i=l---N 



< 



N 



The propagation of chaos property (iii) stems from the almost sure convergence 
of (V^^[t),. . . ,V^^[t),t < T) towards (y^^ (t), (t), t < T), as a process and 

uniformly for fixed time, and is proved in a similar fashion. □ 



The P equations (1.2 1, which arc P implicit stochastic differential equations. 



describe the asymptotic behavior of the network. However, the characterization and 
simulation of their solutions is a challenge. Fortunately, due to their particular form 
in our setting, these equations can be substantially simplified, thanks to the following 
proposition. 

Proposition 1.2. Let us assume that = (y^)a=i...p is a P-dimensional 
Gaussian random variable. We have: 



The solutions of the P mean field equations (1.2 1 with initial conditions 
are Gaussian processes for all time. 

Let fj,{t) — (^a(i))a=i...p denote the mean vector of the process {Va)a=i...p 
andv{t) = {va(t))a=i...p its variance. Let also fi3{x,y) denote the expectation 
of Sp{U) for U a Gaussian random variable of mean x and variance y. We 
have: 



Va{t) 



a = l. 



with initial condition iJ,a{0) = E [V^] and Vq{0) — E [(^^^ — ^0(0))^ 
equation (1.5), the dot denotes the differential with respect to time. 

Proof. 



..P 
(1.5) 
1. In 



The unique solution of the mean field equations (1.2) starting from a 
square integrable initial condition y° measurable with respect to J- can be written in 
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the form: 



J2Jc.i3E[SpiVpis))])ds 

13=1 



(1.6) 



This fixed point equation is shown in 49] to have a unique solution, which is also the 
unique limit of the iteration of the map $ = ($(•)", a = 1, • • • ,P), defined on the 
space of stochastic processes: 



{la 



/3=1 



In other words, for any square-integrable stochastic process X, the unique solution of 
the mean field equations is equal in law to the limit of <i>"(X) as n goes to infinity. 
We also observe that if X° is a Gaussian random variable, $(X) is also Gaussian 
whatever X. We hence conclude that the unique solution of the mean field equation 
is Gaussian, as a limit in law of a sequence of Gaussian processes. 

The law of this solution is hence characterized by its mean and covariance func- 
tions. The formula (1.6) involves the expectation E [S'^(V^(s))], which, because of 
the Gaussian nature of Xp, only depends on fJ.p{s) and vp{s), and is denoted by 
fp{lip{s), vp{s)). Taking the expectation of both sides of the equality ( 1.6 ), we obtain 
the equation satisfied by the mean of the process iJLa{t) — E [V]a(i)]: 



E[F„(0)] 



\J^pE[Sp{Vp{s)]+Io.{s) ds 



Taking the variance of both sides of the equality (1.6), we obtain the following equa- 
tion: 



Va{t) = e (^fa(O)-H j e-^Xlis)dsj, 

and this concludes the proof. □ 
Remark. 

• In order to fully characterize the law of the process V , we need to compute the 
covariance matrix function Cafi{ti,t2) for ti and t2 inH'^*. It is easily shown that: 



tiAt2 '(Tg+Tg) 

e """-s \^{s)\p{s)ds (1.7) 



for t\,t2 G R , hence only depends on the parameters of the system and is in 
particular not coupled to the dynamics. The description of the solution given by 
equations (1.51 is hence sufficient to fully characterize the solution of the mean field 



equations ( 1.2 1 
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The uniformity in time of the propagation of chaos has deep implications in regard 



of equations (1.51. Indeed, we observe that the solution of the mean field equation is 
governed by the mean of the process, the expectation being a deterministic function 
depending on the parameters of the system. The uniformity in particular implies 
that, for i m population a: 

sup |e W^] - < supE \\V: - V:\] < ^= (1.8) 

t>n I L J I t>n L J V Ai 



implying uniform convergence of the empirical mean, as a function of time, towards 

fla{t). 



If V'^ is not Gaussian, the solution of equation (1.2 1 asymptotically converges ex- 
ponentially towards a Gaussian solution. The important uniformity convergence 
property towards the mean field equations ensures that the Gaussian solution is in- 
deed the asymptotic regime of the network, which strengthens the interest of the 



analysis of the differential .system (1.51 



The functions depend on the choice of the sigmoidal transform. A particularly 
interesting choice is the erf sigmoidal function Sa{x) — erflgaX + ja)- In that case 
we are able to express the function in closed form, because of the following lemma: 

Lemma 1.3. In the case where the .sigmoidal transforms are of the form Sa{x) = 



erf(gca; + 7Q,), the functions faili-a^Va) involved in the mean field equations (1.5) with 
a Gaussian initial condition take the simple form: 

/„(M,.)=erf f^^i^y (1.9) 
Proof. The proof is given in Appendix □ 

In summary, we have shown that, provided that the initial conditions of each 
neuron are independent and Gaussian, the large-scale behavior of our linear model 
is governed by a set of ordinary differential equations (theorem |1.1| and proposition 



1.2). This is very interesting since it reduces the study of the solutions to the very 



complex implicit equation (1.2 1 bearing on the law of a process to a much simpler 
setting, ordinary differential equations. As shown below this allows us to understand 
the effects of the system parameters on the solutions. For this reason we assume from 
now on that the initial condition is Gaussian, and focus on the effect of the noise on 
the dynamics. 

2. Noise-induced phenomena. In this section we mathematically and numer- 
ically study the influence of the noise levels Aq, on the dynamics of the neuronal 
populations. Thanks to the uniform convergence of the empirical mean towards the 



mean of the mean field system (equation ( 1.8 )) and the propagation of chaos property 



for the network process, it is relevant to study such phenomena through the thorough 



analysis of the solutions of the mean field equations given by the ODEs (1.5). This is 
what we do in the present section. 



As observed in equation (1.5), in the case of a Gaussian initial condition, the 



equation of the variance v is decoupled from the equation on the mean /i in (1.5). 
The variance satisfies an autonomous equation: 

2 

Va = Va + >?a{t)- 

Ta 
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which is easily integrated as: 

Vctit) is therefore independent of the mean /x. This implies that the equations on jj 
are a set of non-autonomous ordinary differential equations. 

These ordinary differential equations arc similar to those of a single neuron. They 
differ in that the terms in the sigmoidal functions depend on the external noise levels 
Xa{t)- They read: 



p 



Aa = h Ja^eif 



l+gle-^i/-^{vp{Q)+jle^^/-^\l{s)ds} 




hence the slope gp and the threshold 7^ are scaled by a time- varying coefficient which 
is always smaller than one. 

We now focus on the stationary solutions when the noise parameter A does not 
depend upon time. In that case, the variance is equal to: 

Va{t) = Tc,>?j2 + e-^K(O) - T„A^/2), 

and converges exponentially fast towards the constant value TaA^/2. Asymptotic 
regimes of the mean field equations are therefore Gaussian random variables with 
constant standard deviation. Their mean is solution of the equation: 



I,--- ,P 



In other words, the presence of noise has the effect of modifying the slope <?« and 
the threshold 7^ of the sigmoidal erf function, but the type of the equations is the 
same as that of the equation of each individual neuron or cortical column, it is a rate 
equation. 

We observe that the larger the noise amplitude A, the smaller the slope of the 
sigmoidal transform. Noise has the effect of smoothing the sigmoidal transform. This 
will have a strong influence on the bifurcations of the solutions to the mean field 
equations and hence on the behaviors of the system. We demonstrate these effects for 
two simple choices of parameters in one- and two-population networks. 

2.1. The external noise can destroy a pitchfork bifurcation. Let us start 
by considering the case of a one-population network. We drop the index a since no 
confusion is possible. We assume for simplicity that the threshold 7 of the sigmoid is 
mill and that the time constant t is equal to one. By doing so, we do not restrict the 
generality of the study, since r can be eliminated by rescaling the time and 7 can be 
absorbed into / by a simple change of origin for /i. The network equations read: 



dv^ =\-vi + ^Y. {9 vi) + m dt + xdBi z = 1, • 



Noise effects in mean field dynamics 



11 



and we are interested in the limit in law of their solutions as the number of neurons 
N tends to infinity. 

In order to analytically study the effect of the parameter A, we set I = 
In that case, and in the absence of noise, the solution F = is a fixed point of 
the network equations. The following proposition characterizes their solutions in the 
deterministic and stochastic cases. 

Proposition 2.1. In a non- stochastic finite-size network, the null solution is: 

• stable if J < or if J > and g < Qc '■— \/2tt/ J, 

• unstable for J > and g > gc- 

• For J > 0, the system undergoes a pitchfork bifurcation at g — ga- 
in the mean field limit of the same stochastic network, the pitchfork bifurcation occurs 
for a new value of g ^ g* = ^^^^^2 ' > 9c if J > and A < J/^/tt. Furthermore the 
null solution is: 

• stable if: 

— J <0 or 

— J > and A > J/\/t^ (large noise case) or 

— J > 0, X < J/^/tt and g < g* , 

• unstable for J > 0, X < J/\^ and g > g* , and 

• the system undergoes a pitchfork bifurcation at g = g* when J > and 
X < J/v^. 

This proposition is a bit surprising at first sight. Indeed, it says that noise can 
stabilize a fixed point which is unstable in the same non-stochastic system. Even 
more interesting is the fact that if the system is driven by a sufficiently large noisy 
input, the zero solution will always stabilize. It is known, see, e.g., |34| . that noise 
can stabilize the fixed points of a deterministic system of dimension greater than or 
equal to 2. The present observation extends these results to a one-dimensional case, 
in a more complicated setting because of the particular, non-standard, form of the 
mean field equations. Also note that this proposition provides a precise quantification 
of the value of the parameter that destabilizes the fixed point. This is a stochastic 
bifurcation of the mean field equation (a P-bifurcation -P for phenomenological- in 
the sense of P). This estimation will be used as a yardstick for the evaluation of the 
behavior of the solutions to the network equations in section [3j 

Proof. We start by studying the finite-size deterministic system. In the absence 
of noise, it is obvious because of our assumptions that the solution = for all 
i G {1, . . . , TV} is a fixed point of the network equations. At this point, the Jacobian 
matrix reads —Idj\f + j^ Ijv, where Idj^i is the NxN identity matrix and l^v is the 
N X N matrix with all elements equal to one. The matrix Ijv is diagonalizable, all its 
eigenvalues are equal to zero except one which is equal to N. Hence, all eigenvalues of 
the Jacobian matrix are equal to —1, except one which is equal to The solution 

where all are equal to zero in the deterministic system is therefore stable if and 
only a g J < \/2tt. The eigenvalue corresponding to the destabilization corresponds 
to the eigenvector 1 whose components are all equal to 1. Interestingly, this vector 
does not depend on the parameters, and therefore it is easy to check that at the point 
g — gc the system loses stability through a pitchfork bifurcation. Indeed, because of 
the symmetry of the erf function, it appears that the second derivative of the vector 
field projected on this vector vanishes, while the third derivative does not (it is equal 
to-(l+5')). 

Considering now the stochastic mean field limit, the stationary mean firing rate 
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in that case is solution of the equation: 

^l = -/i + Jcrf \ , I + / 

Here again, the null firing rate point /i = is a fixed point of the mean field equations, 
and it is stable if and only if — 1 + J . ^ < 0. The remaining of the 

proposition readily follows from the fact that the stability changes at g = g* where 
J , 1. □ 

727r(l+r^AV2) 

Note that the results in this proposition only depend on A and its effect on the 
slope of the sigmoid. It is a general phenomenon that goes beyond the example in this 
section: increasing A decreases the slope of the sigmoidal transform and the threshold. 
In section |3] we will see that this phenomenon can be observed at the network level, 
and a good agreement will be found between the finite-size network behavior and the 
predictions obtained from the mean field limit. 

We now turn to an example in a two-dimensional network, where the presence of 
oscillations will be modulated by the noise levels. 

2.2. The external noise can destroy oscillations. The same phenomenon 
of nonlinear interaction between the noise intensity and the sigmoid function can 
lead, in higher dimensions, to more complex phenomena such as the disappearance or 
appearance of oscillations. In order to study phenomena of this type, we instantiate a 
simple two-populations network model in which, similarly to the one-dimensional case, 
all the calculations can be performed analytically. The network we consider consists 
of an excitatory population, labeled 1, and an inhibitory population, labeled 2. Both 
populations are composed of the same number N/2 of neurons (A^ is assumed in all 
the subsection to be even), and have the same parameters ti = T2 = r, gi — g2 ^ g 
and Ai A2 ~ A. We choose for simplicity the following connectivity matrix: 

and we assume that the inputs are set to /i = and I2 = — J- The zero solution 
where all neurons have a zero voltage is a fixed point of the equations whatever the 
number of neurons N in each population. We have the following result: 

Proposition 2.2. In the deterministic finite-size network, the null solution is: 

• stable if J < or if J > and g < gc '■= V^tt/ J, 

• unstable for J > and g > gc and the solutions are oscillating on a periodic 
orbit. 

• For J > the system undergoes a supercritical Hopf bifurcation at g = gc- 
In the mean field limit of the same stochastic network, the Hopf bifurcation occurs for 
a new value of the slope parameter g — g* = ^jf^^^xi 9c- Furthermore the null 
solution is: 

• stable if: 

- J <0 or 

- J > and A > J/^/tt ("large noise case") or 

- J > 0, A < J/y/n and g < g* , 

• unstable for J > 0, A < J/\/t^ and g > g* , and the system features a stable 
periodic orbit. 
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• The system undergoes a supercritical Hopf bifurcation at g — g* when J > 
and A < J/y^. 

Note that proposition |2 . 2| is quite similar to proposition |2.1[ the qualitative differ- 
ence being that the system is oscillating. The proof is closely related and is presented 
in less details. 

Proof. In the deterministic network model, the Jacobian matrix at the null equi- 
librium can be written as 



A = -Mn + M<S)1n/2 



9_ 

where (g) denotes the Kronecker product (see e.g. [321 [H]), i-S- the Jacobian matrix is 
built from N/2 blocs of size 2x2 and each of these blocks is a copy of M. The 
eigenvalues of a Kronecker product of two matrices are all possible pairwise products 
of the eigenvalues of the matrices. Since the eigenvalues of M are equal to ^ (1 ± i) 
where — —1, and as noted previously, the eigenvalues of lAr/2 are with multiplicity 
N/2 — 1 and N/2 with multiplicity 1, we conclude that the Jacobian matrix A has 
N — 2 eigenvalues equal to —1, and two eigenvalues equal to —1+g J/V2n{l ± i). 
The null equilibrium in this deterministic system is therefore stable if and only if the 
real parts of all eigenvalues are smaller than 0, viz. gj < \/2i. Therefore, for a fixed 
J, the system has a bifurcation at g^ = \/2j^ / J . The analysis of the eigenvectors 
allows to check the genericity and transversality conditions of the Hopf bifurcation 



(see e.g. [29]) in a very similar fashion to the proof of proposition 2.1 

In the mean field model, the same analysis applies and, as in the one-dimensional 
case, the bifurcation point is shifted to g* when this value is well-defined, which 



concludes the proof of the proposition 2.2 □ 



We have therefore shown that noise can destroy the oscillations of the network. 



The results of propositions 2.1 and 2.2 are summarized in figure 2.1 

An even more interesting phenomenon is that noise can also produce regular 
cycles in the mean part of the solution of the mean field equations, for parameters 
such that the deterministic system presents a stable equilibrium. This is the subject 
of the following section. 

2.3. The external noise can induce oscillations. In order to uncover further 
effects of the noise on the dynamics, we now turn to the numerical study of a slightly 
more complex two-populations network. As in the previous sections, the sigmoids 5*1 
are 5*2 are assumed to be erf functions with the same slope g. The related nonlinear 
functions f\ and /2 involved in the mean field equations are then equal, given by 



direct application of lemma 1.3 and noted /. Both populations receive an additive 



noise with the same intensity A, independent of time. The connectivity matrix is: 

(2.) 



where j accounts for the total connectivity strength of the network. 

This system was studied in a different context in 47 where it displayed an in- 
teresting, non-trivial dynamics. We further assume that the input currents 1\ and li 
do not depend on time and that the initial condition on the variance is the same for 
each population with the effect that the variances of the two populations are equal. 
We denote this variance by v(t). 
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\^ J 

Figure 2.1. Summary of the results in propositions \2. i\ and \2^ see text. The additive noise 
parameter A smoothly modifies the pitchfork or the Hopf bifurcation curve in the {g, J) plane. For A 
large enough, the null solution of the mean field equation is always stabilized whatever g.MF: mean 
field limit, Deterministic: finite-size deterministic network. The blue (respectively red) curve is one 
branch of the hyperbola of equation gj = \/2tt (respectively g\J — ttA^ = \/27rJ. 



The mean field equations in that case read: 

V = -2i^ + A2 

V. r 

The input currents /i and I2 are key to understanding the behavior of the solu- 
tions. Their interaction with the additive noise A profoundly modifies this behavior, 
as we show now. The road to understanding this interaction is to study the bifurca- 
tions of codimension one and twc[^ of the solutions to the mean field equations with 
respect to the parameters Ii, I2, J and A. 

2.3.1. Noise and input. We start by investigating the system for a fixed value 
of the total connectivity strength j — I and study the influence of the noise on the 
behaviors as an input parameter is varied. More precisely, we want to define classes 
of behaviors corresponding to different bifurcation diagrams, as the deterministic ex- 
ternal input parameter Ii is left free while the other parameter I2 is kept constant 
and equal to —3. We proceed with the results of the numerical analysis of the bifur- 
cations of the system, using XPPAut for codimension one diagrams and MatCont 
package [51] 120] for codimension two diagrams. 

Qualitative differences in the dynamics of the neuronal population in response to 
different inputs appear when varying the noise level A. Let us consider the bifurcation 
diagrams of the mean field equations as Ii is varied, for different values of A (the 
qualitative results turn out to change smoothly when I2 is also allowed to vary). 



■tit turns out that there are no codimension three bifurcations. It is therefore sufficient to study 
the bifurcations of codimension one and two. 
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Figure 2.2. Codimension two bifurcation diagram and zooms for the mean field equations as 
II and A are varied. We observe 6 qualitatively different zones (A) through (E), and for legibility 
merged two zones (CI) and (C2). Blue: saddle-node bifurcations, pink: H op f bifurcations, green: 
saddle homoclinic bifurcations, BT: Bogdanov Takens bifurcation, CP: cusp. Individual behaviors 
in each zone are summarized in appendix[B\ 



The corresponding codimension two diagram is shown in Figure |2.2[ We observe 
that the system features three codimension two bifurcations: two cusps (CP) and 
one Bogdanov- Takens bifurcation (BT). In addition to these local bifurcations, we 
observe that the Hopf bifurcation manifold (shown in pink in Figure 2.2 1 has a turning 
point, i.e. changes monotony as a function of A. Similarly, the saddle-homoclinic 
bifurcation curve (shown in green in Figure 2.2 1 has a turning point. Numerical 
values for these turning points are provided in appendix [B] Further analysis of this 
bifurcation diagram shows that there exist six different noise levels that correspond 
to qualitatively different bifurcation diagrams as Ii is varied. Each A-parameter zone 
is labeled in Figure [2?2] with the letters (A) through (E). Two zones, (CI) and (C2), 
are merged into one zone (C) for legibility of the figure. A precise description of the 
behavior of the system in each of these zones is provided in appendix [B] 

We observe that for A large enough no cycles exist for all inputs Ii (zones (D-E)). 
When A is small enough (zones (A) and (B)), cycles always exist: they appear though 
a supercritical Hopf bifurcation and disappear through a homoclinic bifurcation. Zone 
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Figure 2.3. Creation and destruction of cycles by noise as shown in the codimension 1 bifurca- 
tion diagram as a function of X for I\=0. It features three bifurcations: a saddle-node bifurcation 
(LP), a Hopf bifurcation (H) and a saddle homoclinic bifurcation (green circles). We observe three 
main different noise regimes: a high-state equilibrium regime, a periodic regime and a low-state 
equilibrium regime. A small interval of values of A corresponds to the co-existence of cycles and a 
fixed point close to the saddle-homoclinic orbit. 



(C) is a zone where the system is very sensitive to the variations of A. In this zone, 
the Hopf and the homochnic bifurcation curves (arising from the Bogdanov Takens 
bifurcation) have a turning point. Therefore, in this very small parameter zone cor- 
responding to A € [2.935,2.965], the system can either have two Hopf bifurcations 
and two saddle homoclinic bifurcations (CI), i.e. two disconnected branches of cycles 
that appear through supercritical Hopf bifurcations and disappear through homoclinic 
bifurcations, or two supercritical Hopf bifurcations whose family of limit cycles are 
identical (C2). The different behaviors are displayed and described in appendix [b| 

We can also use this diagram to investigate the behavior of the system for a fixed 
value of 1 1 when the noise intensity A is varied. In a similar fashion to what was 
done for fixed A, we observe different ranges of values of Ii where the qualitative 
behavior of the system as a function of A does not vary. Let us for instance fix Ii = 0. 
The corresponding codimension one bifurcation diagram is shown in Figure |2.3[ We 
observe that for small noise intensities, the system features a stable fixed point. This 
fixed point looses stability through a saddle-node bifurcation. A saddle-homoclinic 
bifurcation occurs on the branch of unstable fixed points, and cycles of large amplitude 
and small frequency appear. This phenomenon can be compared to epileptic spikes: 
global oscillations of large amplitude and small frequency suddenly appearing in a 
population of neurons (see PS^). This comparison appears to be also relevant from 
the microscopic viewpoint: in the network simulations of section|3j we observe that all 
neurons suddenly synchronize in these parameter regions, a characteristic phenomenon 
occurring in epileptic seizures. The amplitude of these cycles progressively decreases 
and their frequency progressively increases as the noise intensity is increased, and 
they eventually disappear through a supercritical Hopf bifurcation. 
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2.3.2. Noise and total connectivity strength. Similarly to the previous case, 
we now consider the dynamics of the mean field system as a function of the connectiv- 
ity strength j and the noise, for a fixed value of the currents Ii = and I2 = —3. The 
total connectivity strength j is a particularly interesting parameter for applications, 
for it can account for such phenomena as the well established fact that functional 
connectivity is increased in epilepsy (see e.g. [8]). This influence was studied in an- 
other context by Pham and collaborators in ^T]. In this contribution, the authors 
study randomly or fully connected networks of spiking neurons. They analyze the 
probability distribution of spike sequences and reduce this analysis to the study of 
the properties of a certain map under an independence assumption and in the limit 
where the number of neurons is infinite (which makes the independence assumption 
particularly relevant). They show that noise can trigger oscillations for certain values 
of the total connectivity parameter j in a one-population case. They end up with a 
partition of the parameter space in different zones where the system either shows a 
single "high" fixed point, a "low" fixed point or oscillations, or multistability between 
these different attractors. 



In order to compare ours to these previous results, we adopt the same presen- 
tation in this section. The codimension two bifurcation diagram with respect to j 
and A contains manifolds of saddle-node, Hopf and saddle-homoclinic bifurcations, 
and shows no codimension two bifurcation. The parameter space can be partitioned 
into five qualitatively different zones where the system either shows a "low" fixed 
point (L), a "high" fixed point (H), oscillations (O), or multistability between these 
attractors (see Figure 2.4 1. The diagram is quite similar to the one found in As 
in the case discussed in section |2.3.H the analysis of this diagram results in defining 
three regions of values of the noise parameter A in which the bifurcation diagrams 
as a function of j are qualitatively identical. These three regions make less sense 
here as the total connectivity parameter is not easily controllable experimentally, in 
contrast with the external input parameter Ii. The resulting diagram has a number 
of similarities and dissimilarities with the results of Pham and collaborators [411 Sec- 
tion V and Fig. 11]. Indeed, similarly to the spiking case, we observe that for small 
values of the connectivity strength j the fully connected system is characterized by 
a low fixed point, whereas large values of the connectivity parameter chiefly corre- 
spond to the presence of a high fixed point, and the overall structure is comparable. 
One of the main differences is that in our case, a two-population network, the system 
does present oscillations as the only attractor for some parameter values. Another 
interesting difference is that for small values of the noise parameter, the spiking sys- 
tem studied in |41j never presents oscillations, whereas the diagram 2.4 exhibits an 
oscillatory (O) and a bistable (OH) regions for arbitrarily small values of the noise 
parameter. This is due to the fact that their system consists of a unique neural popu- 
lation with purely excitatory interactions, whereas our network models the interaction 
of an excitatory and an inhibitory populations, and can feature deterministic oscil- 
lations. We eventually note that in one-population linear firing rate models with 
additive noise and constant external forcing, no oscillatory activity is possible in the 
mean field limit, since its dynamics can be reduced to a one-dimensional autonomous 
dynamical system. Smooth nonlinearities in the intrinsic dynamics or discontinuities 
such as the presence of a spiking threshold in |41j makes the dynamics of the mean 
field equations more complex, in particular prevents reduction to a one-dimensional 
autonomous system governing the mean of the solution. Such intricacies may be the 
source of oscillations in one-population systems. 
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Figure 2.4. Codimension two bifurcation diagram of the mean field system with respect to the 
noise intensity A and the total connectivity j. The diagram is partitioned into 5 different zones 
depending on the number and type of stable attractors: L: Low fixed point, H: High fixed point, O: 
Oscillations. Transitions between two zones are characterized by codimension one bifurcations: Blue 
line: saddle-node bifurcation, red: Hopf bifurcation, green: saddle homoclinic bifurcation. 



We therefore conclude from the analysis of these simulations that noise does 
not only destroy structures and regularity, it can also generate oscillations. These 
noise-induced oscillations are very interesting from the biological viewpoint. Indeed, 
oscillations are essential for the brain function. The link between oscillations and noise 
level is therefore a very relevant piece of information, that strengthens interpretations 
of the functional role of the noise. We develop further these topics in the discussion 
section. 

3. Back to the network dynamics. Thus far, we studied the dynamics of the 
mean field equations representing regimes of the network dynamics in the limit where 
the number of neurons is infinite. We now compare the regimes identified in this 
analysis with simulations of the finite-size stochastic network. We are particularly 
looking for potential finite-size effects, namely qualitative differences between the 
solutions to the network and the mean field equations. This will provide us with 
information about the accuracy of an approximation of the network dynamics by the 
mean field model, as function of the size of the network. 



3.1. Numerical Simulations. Numerical simulations of the network stochastic 



differential equations (1.1) are performed using the usual Euler-Maruyama algorithm 
(see e.g. [351 [M]) with fixed time step (less than 0.01) over an interval [0, T]. In order 
to observe oscillations, we choose T between 50 and 70. The simulations are performed 
with Matlab(R), using a vectorized implementation that has the advantage to be very 
efficient even for large networks. The computation time stays below Is for networks 
up to 2 000 neurons, and appears to grow linearly with the size of the network once 
the cache memory saturates (see Figure 3.1 1. For instance, for T = 20, dt = 0.01, 



the simulation of a 2 000 neurons network takes 0.89s, and for 525 000 neurons, 600s 
on a HP Z800 with 8 Intel Xeon CPU E5520 @ 2.27 GHz 17.4 Go RAM. The main 
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Figure 3.1. Computation time for the simulation of the stochastic network in logarithmic scale 
as a function of the network size. 



limitation preventing the simulation of very large networks is the amount of memory 
required for the storage of the trajectories of all neurons. 

An important property arising from theorem is that asymptotically, neurons 
behave independently and have the same probability distribution. In our numerical 
simulations, we will make use of this asymptotic independence and, in order to eval- 
uate an empirical mean of the process related to a given neuron in population a, will 
compute both the empirical mean over all neurons in that population and a mean 
over different independent realization of the process. This method allows to reduce 
sensitively the number of independent simulations in order to obtain a given precision 
in the empirical mean evaluation. 

3.2. A one population case. We start by addressing the case discussed in 
section 



2.1 where we showed analytically that the loss of stability of the null fixed 
point as the slope of the sigmoid was varied depended on the noise parameter A. We 
now investigate numerically the stability of the fixed point of the network equations. 
In order to check for the presence of a pitchfork bifurcation, we compute, for each 
value of the noise and for each value of the slope of the sigmoid, an estimated value 
of the mean of the membrane potential. This estimate is calculated by averaging 
out over 500 independent realizations the empirical mean of the membrane potentials 
of all neurons in the network at the final time. We display the average value then 
compare these simulations with those of the mean field equations stopped at the same 
time as the network. We observe that both are very similar and show some differences 
with the bifurcation diagram that corresponds to the asymptotic regimes. 

The results of the simulations, where we have also varied N, are shown in Fig- 
ure 3.2 and reveal two interesting features. First, because we simulate over a finite 
time, we tend to smooth the pitchfork bifurcation: this is perceptible for both the 
network and the mean field equations. Second, we observe that the loss of stability of 
the zero fixed point arises at the value of A predicted by the analysis of the mean field 
equations for networks as small as 50 neurons. The value reached by the simulations 
of the network is very close to that related to the mean field equation as soon as N 
becomes greater than 250. 

3.3. Two populations case and oscillations. We now investigate the case 



shown in Figure 2.3 where cycles are created (through homoclinic bifurcation) or 
destroyed (through Hopf bifurcation) as the additive noise intensity parameter A is 
increased. 



Looking at Figure 2.3 we observe that for A £ [1.12, 1.33], stable periodic orbits 
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Figure 3.2. Comparison of the pitchfork bifurcations with respect to the slope parameter g for 
the network and the mean field equations, T = 40, dt = 0.001, number of sample paths: 100, initial 
condition = 0.5 (hence we only see the positive part of the pitchfork, symmetrical solutions 
are found for negative initial conditions, and are not plotted for legibility), (a): 50 000 neurons. 
Continuous curves correspond to network simulations, dashed curves to mean field simulations. 
When A increases, as predicted by proposition \2.1\ we observe that the value of the parameter g 
related to the pitchfork bifurcation increases as well, until the pitchfork disappears: red: A = 0, blue: 
A = 0.4, green: A = 0.8 > \/2n/J ~ 0.56. (b): A = 0.4. The solution to the mean field equation 
undergoes a pitchfork bifurcation at g = 3.55. Large dotted red: theoretical pitchfork bifurcation. 
Large black: endpoint of mean field simulation at time T = 40. The other colored curves show the 
results of the network simulation for different values of the size of the network N . The solution, 
which loses stability, is displayed in thin dashed black. We see that as N increases, the mean field 
equation describes accurately the network activity. For N > 50 (red, green, dotted blue and dotted 
cyan curves) the bifurcation diagram is quite close to the one predicted by the mean field analysis. 



coexist with stable fixed points in the mean field system. For smaller values of A, the 
mean field system features a unique stable fixed point, while for A € [1.33,1.97], it 
features a unique stable limit cycle, and for A > 1.97, the dynamics is reduced to a 
unique attractive fixed point. Numerical simulations confirm this analysis. Let us for 
instance illustrate the fact that the network features the bistable regime, the most 
complex phenomenon. Figure |3.3| shows simulations of a network composed of 5 000 
neurons in each population (time step dt ~ 5 ■ 10"'^, total time T = 50). Depending 
on the mean and on the standard deviation of the initial condition, we observe that 
the network either converges to the mean field fixed point or the periodic orbit. Both 
mean field equations show very close behaviors. 

In the fixed point regime corresponding to small values of A we observe that the 
membrane potential of every neuron randomly varies around the value corresponding 
to the fixed points of the mean field equation (see Figure 3.4 cases (a) and (b))), with 
a standard deviation that converges toward the constant value A^/2 as predicted by 
the mean field equations. The empirical mean and standard deviation of the voltages 
in the network show a very good agreement with the related mean field variables. For 
larger values of A corresponding to the oscillatory regime (Figure 3.4 cases (c) and 
(d)), all neurons oscillate in phase. These synchronized oscillations yield a coherent 
global oscillation of the network activity. The statistics of the network are again 
in good agreement with the mean field solution. The standard deviation converges 
towards the constant solution of the mean field equation. This is visible at the level 
of individual trajectories, that shape a "tube" of solutions around the periodic mean 
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(a) A = 1.2. Oscillatory regime. Statistics (b) A = 1.2. Fixed-point regime Statistics 
of the network compared to the mean field. of the network compared to the mean field. 



Figure 3.3. Featuring bistability. In both cases A = 1.2. The initial conditions for the mean 
field equation are chosen in agreement with the initial conditions of the network. The initial value 
of the membrane potential of each individual neuron in the network is drawn independently from 
a Gaussian distribution of variance 1 whose mean depends on the population: (a) mean 0.5. (b) 
mean 4. Cyan (resp. magenta) curves: value of the mean variable of the mean field solution for 
population 1 (resp. 2). Dashed blue (resp. red) curves: empirical mean of population 1 (resp. 2). 
Yellow: value of the variance of the mean field solution. Dashed black (resp. green): empirical 
variance of population 1 (resp. 2). 



field solution, whose size increases with A. The empirical means accurately match the 
regular oscillations of the solution of the mean field equation. A progressive phase shift 
is observed, likely to be related with the time step dt involved in the simulation. Note 
that the phase does not depend on the realization. Indeed, according to theorem 
the solution of the mean field equations only depends on the mean and the standard 
deviation of the Gaussian initial condition, which therefore governs the phase of the 
oscillations on the limit cycle (see Figure 3.5). 

In the fixed point regime related to large values of A, very noisy trajectories are 



obtained because of the levels of noise involved (see Figure 3.4 cases (e) and (f)) 



Though the individual neurons show very fluctuating trajectories, the empirical mean 
averaged out over all neurons in the network fits closely the mean field fixed point 
solution. 

Eventually, we study the switching between a fixed-point regime and an oscilla- 
tory regime by extensively simulating the 10 000 neurons network for different values 
of A and computing the Fourier transform of the empirical mean (see Figure 3.6). 
The three-dimensional plots show that the appearance and disappearance of oscilla- 
tions occur for the same values of the parameter A as in the mean field limit, and 
the route to oscillations is similar: at the homoclinic bifurcation in the mean field 
system, arbitrarily small frequencies are present, this is also the case for the finite- 
size network. At the value of A related to the Hopf bifurcation, the system suddenly 
switches from a non-zero frequency to a zero frequency in a form that is very similar 
to the network case. Therefore we conclude that the mean field equations accurately 
reproduce the network dynamics for networks as small as 10000 neurons, and hence 
provide a good model, simple to study, for networks of the scale of typical cortical 
columns. As a side remark, we note that at a homoclinic bifurcation of the mean field 
system, very small frequencies appear and a precise description of the spectrum of the 
network activity would require very large simulation times to uncover precisely the 
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(a) A = 0.6. Fixed-point regime. Individ- (b) A = 0.6. Fixed-point regime Empirical 
ual trajectories vs mean field. network statistics vs mean field. 
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(c) A = 1.2. Oscillatory regime. Individ- (d) A = 1.2. Oscillatory regime Empirical 
ual trajectories vs mean field. network statistics vs mean field. 
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(e) A = 2.5. Noisy fixed point regime. (f) A = 2.5. Noisy fixed point regime 

_i 4- — : — 4- — : field. Empirical network statistics vs mean field. 



Individual trajectories 



vs mean 



Figure 3.4. Solution of the network dynamics for different values of the noise parameter A 
compared to the mean field solution. Simulations are run for 10 000 neurons, 5 000 in each popula- 
tion, (a), (c), (e): 40 individual trajectories of the membrane potentials of AO neurons arbitrarily 
chosen in the network (20 in each population) compared to the solution of the mean field equations. 
Blue: population 1 (excitatory) . Red: population 2 (inhibitory). Cyan (resp. magenta): mean of 
the mean field solution for population 1 (resp. 2). (b), (d), (f): Empirical statistics of the network 
compared to the mean field. Cyan (resp. magenta): mean of the mean field solution for population 
1 (resp. 2). Yellow: variance of the mean field solution. Dashed blue (resp. red): empirical mean 
of population 1 (resp. 2). Dashed black (resp. green): empirical variance of population 1 (resp. 2). 
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Figure 3.5. Different realizations of the stochastic network dynamics: the membrane potentials 
of 5 neurons among 5 000 of population 1 are plotted for 12 different realizations represented in 
different colors. All neurons oscillate in phase, and this phase does not depend on the realization. 



spectrum at this point, even more so since the large standard deviation of the process 
disturbs the synchronization. In a forthcoming study focusing on mean field equations 
arising in a similar system including dynamically fluctuating synaptic coefficients, an 
interesting additional phenomenon will appear: in that case, the standard deviation 
variable will be coupled to the mean, and this coupling will result in sharpening the 
synchronization. 

We conclude this section by discussing heuristic arguments explaining the ob- 
served regular oscillations. Let us start by stating that this phenomenon is a pure 
collective effect: indeed, two-neurons networks (one per population) do not present 
such regular oscillations as noise is varied. We observe that individual trajectories 
of the membrane potential of a 2-neurons networks for small noise levels stay close 
to the deterministic fixed point. However, when noise is increased, the system starts 
making large excursions with a typical shape resembling the cycle observed in the 
mean field limit, and these excursions occur randomly. Such excursions are typical 
of the presence of a homoclinic deterministic trajectory: when perturbed, the system 
catches the homoclinic orbit responsible for such large excursions. The codimension 
one bifurcation diagram of the 2-neurons system indeed illustrates the presence of a 
homoclinic orbit as a function of /i (see diagram 2.2 and Figure [bj] (A)|^ Noise 



can be heuristically seen as perturbing the deterministic value of Ii . For sufficiently 
small values of the noise parameter, the probability of Ii to visit regions correspond- 
ing to the presence of a cycle is small. But as the noise amplitude is increased, this 
probability becomes non- negligible and individual trajectories will randomly follow 
the stable cycle. Such excursions produce large input to the other neurons which will 
either be inhibited or excited synchronously at this time, a phenomenon that may 
trigger synchronized oscillations if the coupling is strong enough and the proportion 
of neurons involved in a possible excursion large enough. If the noise parameter is too 
large, the limit cycle structure will be destroyed. 



§ Indeed, the mean field equations with A = are precisely the equations of a two-neurons network 
since in that case f{pL, = S{^i). 
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(c) 



Figure 3.6. Squared moduli of the Fourier transforms of (a) the empirical mean for simulations 
of the network and (b ) of the mean variable of the solution to the mean field equations as functions 
of the frequency (Hz) and the noise parameter A. We observe that oscillations appear in the network 
for the same value of A as in the mean field equations (Figure \2.0^ , first through what appears to 
be a homoclinic bifurcation (arbitrary small frequencies) and also disappear for the same value of \ 
through what seems to be a Hopf bifurcation (discontinuity in the power spectrum), (c) Magnitude of 
the difference between the two diagrams: we note that the frequency distribution reaches its maxima 
for these same values of A, and the main differences are observed, as expected, around the putative 
homoclinic bifurcation point. 



Another way to understand this phenomenon consists in considering the phase 



plane dynamics of the two-neurons network with no noise (see Figure 3.7). The 
system presents three fixed points, one attractive, one repulsive, and a saddle. The 
unstable manifold of the saddle fixed point connects with the stable fixed point in 
an heteroclinic orbit. The stable manifold of the saddle fixed point is a separatrix 
between trajectories that make small excursions around the stable fixed point, and 
those related to large excursions close to the heteroclinic orbit. As noise is increased, 
the probability distribution of each individual neuron, centered around the stable fixed 
point, will grow larger until it crosses the separatrix with a non-negligible probability, 
resulting in the system randomly displaying large excursions around the heteroclinic 
cycle. The fact that a homoclinic path to oscillations is found in the mean field 
limit can be accounted for by these observations, considering the fact that crossing 
the separatrix, when noise is of small amplitude, can take an arbitrary long time. 
The rhythmicity of the oscillations we found and the synchronization are related to 
the coupling in a complex interplay with the probability of large excursions. These 
heuristic arguments require a thorough mathematical analysis that is outside the scope 
of the present paper. 

4. Discussion. In this article, we have been interested in the large-scale be- 
havior of networks of firing rate neuron models. Using a probabilistic approach, we 





Figure 3.7. Trajectories in the phase plane for different values of X superimposed on the phase 
diagram. Red curve: fii-nullcline, Green curve: fi2-nullcline, Orange cycle: unstable manifold of 
the saddle fixed point (heteroclinic orbit) and Cyan curve: stable manifold of the saddle fixed point 
(note that it is almost superposed with part of the /ii-nullcline), constituting the separatrix between 
those orbits that directly return to the stable fixed point and those following the heteroclinic cycle. 
Black: noisy trajectories. Upper left: X = 0.2,- no excursion, corresponds to the fixed point regime. 
Upper right: A = 1.- rare excursions do occur, corresponding to the bistable regime. Bottom left: 
X = 1.6; excursions are frequent but occur irregularly (corresponding to the oscillatory regime). 
Bottom right: A = 5; the heteroclinic cycle structure is lost, corresponding to the fixed point regime. 



addressed the question of the behavior of neurons in the network as its size tends to 
infinity. In that hmit, we showed that all neurons behaved independently and satisfied 
a mean field equation whose solutions are Gaussian processes such that their mean 
and variance satisfy a closed set of nonlinear ordinary differential equations. Uniform 
convergence properties were obtained. 

We started by studying the solutions of the mean field equations, in particular 
their dependence with respect to the noise parameter using tools from dynamical 
systems theory. We showed that the noise had non-trivial effects on the dynamics 
of the network, such as stabilizing fixed points, inducing or canceling oscillations. A 
codimension two bifurcation diagram was obtained when simultaneously varying an 
input parameter and the noise intensity, as well as simultaneously varying a total input 
connectivity parameter and the noise intensity. The analysis of these diagrams yielded 
several qualitatively distinct codimension one bifurcation diagrams for difi'erent ranges 
of noise intensity. Noise therefore clearly induces transitions in the global behavior of 
the network, structuring its Gaussian activity by inducing smooth oscillations of its 
mean. 

These classes of behaviors were then compared to simulations of the original finite- 
size networks. We obtained a very good agreement between the simulations of the 
finite-size system and the solution of the mean field equations, for networks as small as 
a few hundreds to few thousands of neurons. Transitions between different qualitative 
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Figure 4.1. Empirical distribution of the values of (V(T))i=i.. jy /<"" N = 1000 (500 neurons 
per population) in each population (blue and green filled distribution) versus theoretical mean field 
distribution. The Kolmogorov-Smirnoff validates the fit of the distributions (see text). 

behaviors of the network matched precisely the related bifurcations of the mean field 
equations, and no qualitative systematic finite-size effects were encountered. More- 
over, it appears that the convergence of the solution to a Gaussian process as well as 
the propagation of chaos property happen for quite small values of N, as illustrated in 
Figure [4?T| This figure represents the distribution of the voltage potential at a fixed 
time r = 40 for N = 500, simulated for 20 sample trajectories. The Kolmogorov- 
Smirnov test validates the Gaussian nature of the solution with a p-value equal to 
7 • 10^"*. In order to test for the independence, we used the Pearson, Kendall and 
Spearman tests of dependence. We obtain the correlation values 0.0439 (p-value 0.33) 
for the first population, 0.0212 (p-value 0.4785) for the second, and 0.0338 (p-value 
0.45) for the cross-correlation between populations, all of them clearly rejecting the 
dependence null hypothesis. This independence has deep implications in the efficiency 
of neural coding, a concept that we will further develop in a forthcoming paper. 

These findings have several implications in neuroscience and other scientific do- 
mains as we now comment. 

4.1. Noise-induced phenomena. We have seen that the presence of noise in 
the system induces different qualitative behaviors. For instance, regular oscillations 
of the mean firing rate, linked with a synchronization of all neurons in the network, 
appear in the system at some precise values of the noise parameter, in particular 
for systems that feature a stable fixed point in a noiseless context. This means that 
noise has a strong structuring effect on the global behavior of a cortical assembly, 
which is rather a counterintuitive phenomenon, since noise is usually chiefly seen as 
altering structured responses. This phenomenon adds to a few recent observations 
made in other settings by various authors. We related our results to the phenomenon 
of coherent oscillations in one-population spiking neural networks in section [2.3. 2[ in 
relationship with the works of Pham and collaborators ^41j . The results obtained by 
Nesse and collaborators [38] also confirm the presence of noise-induced phenomena 
in all-to-all coupled networks of leaky integrate-and-fire neurons with filtered noise 
and slow activity-dependent currents, i.e. noise and spikes are filtered by a second 
order kernel. The mean field equations they derive are based on a separation of time 
scale and ergodicity properties. They show in that case the presence of noise-induced 



Noise effects in mean field dynamics 



27 



burst oscillations similar to the case of [4ll, and a resonance phenomenon. Both 
approaches differ from our analysis in that they are dealing with spiking neurons. 
Pham and collaborators are able to reduce the dynamics a discrete-time network 
of pulse-coupled spike response model neurons, using a Markovian approach, to the 
study of a discrete time dynamical system. Nesse and collaborators derive their mean 
field limit using relevant approximations of the Fokker-Planck equations. The interest 
of our firing-rate model is that it rigorously allows the use of the mathematical theory 
of propagation of chaos to reduce the system to a set of coupled differential equations 
that can be studied using the dynamical systems theory. 

The phenomena observed in our analysis of large-scale neuronal networks differs 
from the phenomena of stochastic resonance or coherence resonance well documented 
in the neurocomputational literature (see e.g. [3^ for a review of the effect of noise 
in excitable systems). These phenomena correspond to the fact that there exists a 
particular level of noise maximizing the regularity of an oscillatory output related 
to periodic forcing (stochastic resonance) or to intrinsic oscillations (coherence reso- 
nance). Such situations are evidenced through the computation of the maximal value 
of the Fourier transform of the output. In our case, as we can see in the Fourier trans- 
form plots, the maximal value of the Fourier transform does not present a clear peak 



as a function of the noise level (see Figure 3.61, hence the system does not exhibit 
resonance. Besides this observation, the regularity of the oscillation can be expected 
to be relatively high for large networks in our framework, since the mean activity is 
asymptotically perfectly periodic. 

Our study possibly opens the way to applications in different scientific domains 
such as economics where global behaviors emerge from the noisy interaction of dif- 
ferent individual agents, or physics, in particular in thermodynamics where global 
behaviors arise from the interaction of a large amount of particles. 

Another important result of ours is the ability to define classes of parameter 
ranges attached to a few generic bifurcation diagrams as functions of the input to a 
population. This property suggests further some reverse-engineering studies allowing 
to infer from measurements of the system responses to different stimuli the level of 
noise it is submitted to. 

4.2. Understanding the functional role of noise. The question of the func- 
tional role of noise in the brain is widely debated today since it clearly affects neuronal 
information processing. A key point is that the presence of noise is not necessarily 
a problem for neurons: as an example, stochastic resonance helps neurons detecting 
and transmitting weak subthreshold signals. Furthermore neuronal networks that 
have evolved in the presence of noise are bound to be more robust and able to explore 
more states, which is an advantage for learning in a dynamic environment. 

The fact that noise can trigger synchronized oscillations at the network level en- 
riches the possible mechanisms leading to rhythmic oscillations in the brain, directly 
relating it to the functional role of oscillations. Rhythmic patterns are ubiquitous in 
the brain and take on different functional roles. Among those, we may cite visual 
feature integration |43j . selective attention, working memory. Abnormal neural syn- 
chronization is present in various brain disorders [50 . Oscillations themselves can 
signal a pathological behavior. For instance, epileptic seizures are characterized by 
the appearance of sudden, collective, slow oscillations of large amplitude, correspond- 
ing at the cell level to a synchronization of neurons, and visible at a macroscopic scale 
through EEG/MEG recordings. This phenomenon is very close to the observation in 
our model that, as noise is slowly increased, the solutions of the mean field equations 
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undergo a saddle-homoclinic bifurcation abruptly yielding large amplitude and small 
frequency oscillations. Such a collective phenomenon resembles epileptic seizures. 
Moreover, in our model and at the microscopic level, these regimes are characterized 
by a sudden precise synchronization of all neurons in the network, consistent with 
what is observed at the cell level in epileptic seizures. 

This mean field model, based on a simple description of neural activity, is therefore 
able to account for complex biologically relevant phenomena, in a mathematically and 
computationally tractable way while including noise effects. This observation suggests 
to use this new model as a cortical mass model and compare it to more established 
cortical column models such as Jansen and Rit's or Wendling and Chauvel's [311 EH 

4.3. Perspectives. Several extensions of the present study with applications in 
neuroscience and in applied mathematics are envisioned. 

The first is to expand our work to include more biologically relevant models and 
to study the behavior of the solutions of the mean field equations in that mathemat- 
ically much more complex setting that would include in particular nonlinear intrinsic 
dynamics and different ionic populations. Another important direction in the de- 
velopment of this work would consist in fitting the microscopic model to biological 
measurements. This would yield a new neural mass model for large scale areas and 
develop studies on the appearance of stochastic seizures and rhythmic activity in re- 
lationship with different parameters of the model, integrating the presence of noise in 
a mathematically and biologically relevant manner. 

From this point of view, the simplicity of the model, specifically the linearity 
of the intrinsic dynamics of each individual neuron, made possible an analytic and 
quantitative analysis of our mean field equations, a particular case of McKean-Vlasov 
equations. The drawback of this simplicity is that it does not represent precisely the 
activity of individual neurons. For realistic individual neuron models the solutions 
to the mean field equations will not, in general, be Gaussian. General approaches 
to study them consist either in studying their properties as random processes, or in 
describing their probability distribution. In the first case, one is led to investigate 
an implicit equation in the space of stochastic processes, and in the second case, one 
is led to study a complex non-local partial differential equation (Fokker-Planck) . In 
both cases one faces a difficult challenge, and the dependency of the solutions with 
respect to parameters is extremely hard to describe. 

Another interesting improvement of the model would consist in considering that 
the synaptic weights are random. These weights can be randomly drawn in a dis- 
tribution and frozen during the evolution of the network. The propagation of chaos 
would again applies in that case, and in the rather simple case discussed in the present 
manuscript, the solution is a Gaussian process as proved in [27j . However, the mean 
and covariance cannot be described by a set of ordinary differential equations, as in 
here since the covariance depends on the whole previous history of the correlations. 
Alternatively, the weights can be considered as stochastic processes themselves. The 
question of the synaptic noise model chosen and the dynamics of such systems will be 
addressed in a forthcoming paper. 

The present study can be seen as a proof of concept, and it seems reasonable 
to extrapolate that such noise-induced transitions do occur as well for the solutions 
to the mean field equations of these more complex and more biologically plausible 
systems. 
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Appendix A. Proof of lemma |1.3[ 

In this appendix we prove lemma [T3| stating that in the case where the sigmoidal 
transforms are of the form Sa{x) — eTi{gaX + 7a), the functions faiponVa) involved 



form ( 1.9 1 



in the mean field equations (1.2) with a Gaussian initial condition take the simple 



Proof. We have, using the definition of the erf function 



E[S^{X^){t)]^ / erf (5„ (xv/^ai)+/^aW) +7a) 

r-ffo(2;V^„(t)+/^c(t))+7o g-(a-"+y2)/2 



dx 



27r 



dxdy 



This integral is of the form: 



ax+b -(x^+a^)/2 



2tt 



-dxdy 



and therefore, the integration domain has an affine shape as plotted in figure \K7\] In 

y 




Figure A.l. Change of variable for the erf function. 

order to compute this integral, we change variables by a rotation of the axes {x,y) 
and align the affine boundary of our integration domain with our new variables (u, v) 
(see figure A.l). Simple geometric analysis shows that the rotation angle 6 for this 
change of variable is such that tan{9) = a. The new integration domain is in the new 
coordinates given by w < u,„ = bcos{6) 



ax+b -(x^+y'')/2 



2n 



-dxdy 



= erf 



e'^-'+^^'^/'-dudv 



which reads with the parameters of the model: 



faifJ;v) = erf 



Appendix B. Bifurcations Diagram as a function of A. In section [2j we 
observed that six different bifurcation diagrams appear as Ii is varied, depending on 
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Type 


C 


BT 


Horn TP 


H TP 


C 


A 


0.16 


2.934 


2.948 


2.968 


3.74 



Table B.l 

Numerical values of the separation into six X zones for Figure [2. S\ C stands for Cusp, BT: 
Bogdanov-Takens, Horn TP: turning point of the Homoclinic bifurcations curve, H TP: Hopf bifur- 
cation curve turning point. 






^' h ' " ..... . 

Figure B.l. Typical behavior of the system in each zone (A) through (E). (A): A = 0, (B): 
A = 1, (CI): A = 2.945, (C2):\ = 2.955, (D): A = 3, (E): A = 4. Red stars: bifurcations, 
SN: Saddle-Node, H: Hopf, green circle: Saddle- Homoclinic bifurcation. Thick blue line: stable 
fixed point, thin blue line: unstable fixed point, thick pink line: stable cycle. See text for precise 
description. 



the value of A characterizing the additive noise input. For the particular choice of 
parameters chosen in that section, the different zones are segmented for values of A 
given in table |B.l] 

In each of these zones, typical codimension 1 bifurcation diagrams as the input 



Ii is varied are depicted in figure B.l We now describe the behavior of the system 
in each of these zones. 

(A) For very small values of A, the system features four saddle- node bifurcations 
and one supercritical Hopf bifurcation, associated to the presence of stable 
limit cycles that disappear through saddle-homoclinic bifurcation arising from 
the Bogdanov-Takens bifurcation (after the turning point.) In an extremely 
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limited range of parameter, the occurrence of two saddle-node bifurcation 
relates to a bistable regime in that small parameter region. 
(B) In zone (B), the system differs from zone (A) in that the two inferior saddle- 
node bifurcation disappeared through Cusp bifurcation. Globally the same 
behavior are observed, except for the bistable behavior commented above 
(which was not a prominent phenomenon due to the reduced parameter region 
concerned). 

(CI) On the upper branch of saddle nodes, the systems undergoes a Bogdanov- 
Takens bifurcations, yielding the presence in zone (CI) of a supercritical Hopf 
bifurcation and of a saddle-homoclinic bifurcation curve. This BT bifurcation 
accounts for the family of Hopf bifurcations observed in zones (A-B) and for 
the saddle-homoclinic bifurcations, because of the turning points observed 
in the full bifurcation diagrams. In region (CI), two families of limit cycles 
coexist, both arising from supercritical Hopf bifurcation and disappearing 
through saddle-homoclinic bifurcation. 

(C2) Because of the topology of the bifurcation diagram, the turning point of the 
saddle-homoclinic bifurcations curve occurs before the turning point of the 
Hopf bifurcations curve. This difference yields zone (C2) between the two 
turning points. In that zone, we still have two supercritical Hopf bifurcations, 
but no more homoclinic bifurcation. The families of limit cycles corresponding 
to each of the Hopf bifurcations are identical. 

(D) After the turning point of the Hopf bifurcations manifold, we are left with two 

saddle-node bifurcations, hence a pure bistable behavior with no cycle. 

(E) Both saddle-node bifurcation disappear by merging into a cusp bifurcation. After 

this cusp, the system has a trivial behavior, i.e. it features a single attractive 
equilibrium whatever Ii. 
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